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Abstract. The propagation of a superintense laser pulse in an underdense, 
inhomogeneous plasma has been studied numerically by two-dimensional particle-in- 
cell simulations on a time scale extending up to several picoseconds. The effects of the 
ion dynamics following the charge-displacement self-channeling of the laser pulse have 
been addressed. Radial ion acceleration leads to the "breaking" of the plasma channel 
walls, causing an inversion of the radial space-charge field and the filamentation of the 
laser pulse. At later times a number of long-lived, quasi-periodic field structures are 
observed and their dynamics is characterized with high resolution. Inside the plasma 
channel, a pattern of electric and magnetic fields resembling both soliton- and vortex- 
like structures is observed. 
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1. Introduction 

The propagation of superintense laser pulses through low-density plasmas gives rise 
to a variety of nonlinear electromagnetic phenomena [B Ej. As a general issue the 
response of the plasma is nonlinear due both to relativistic effects (hence the definition 
of "relativistic optics" [3]) and to the intense ponderomotive force (i.e., the radiation 
pressure) which strongly modifies the local plasma density. Probably, the example of 
such dynamics which has been mostly investigated is the self-focusing, channeling and 
filamentation of the laser pulse [U El H El [10]. Another prominent effect is the 
generation of coherent structures such as electromagnetic solitons or vortices. Numerical 
simulations (see e.g. [HI [T2], [13]) show that such structures are generated during the 
interaction with the laser pulse on a ultrafast (femtosecond) timescale, but they may 
lead to typical field structures which last for much longer times (e.g. "post-solitons" ) 
[T2] i.e. in the picosecond range, allowing for their experimental observation [H]. On 
such a scale the temporal evolution of such field structures must be studied including the 
effects of the motion of the plasma ions. Stability and evolution of coherent structures 
on the ion time scale has been studied theoretically and numerically in several papers, 
for various regimes and dimensionality [HJ [El EH EE! EE1 CEO HH1 EEH1 [20] . 

In this paper we report a theoretical study of nonlinear effects during and 
after the propagation of a superintense laser pulse in an underdense, longitudinally 
inhomogeneous plasma. The work was motivated by experiments on laser propagation 
in a low-density plasma where the dynamics of self-generated, slowly varying 
electromagnetic fields was investigated using the proton diagnostic technique [2T| . In 
this paper we focus on the simulation results and on their theoretical interpretation, 
while a comparison with the experimental results will be reported elsewhere [22| 123] . 

2. Simulation set-up 

The laser-plasma interaction simulations were performed using a particle-in-cell (PIC) 
code in 2D with Cartesian geometry. Reduction to 2D was dictated by the need to 
address relatively long spatial and temporal scales, close to the experimental ones. 
Moreover (as it will be clear from the discussion of the results) during the interaction 
sharp gradients in the field and current patterns are generated. Thus, a reasonable 
resolution is mandatory to resolve such details, pushing the memory requirements in 3D 
much beyond present-day supercomputing capabilities. Among the set of 2D simulations 
that were performed for the present study, the largest ones employed a 7750 x 2400 grid, 
with spatial resolution Ax = Ay = A/10 (where A is the laser pulse wavelength) and 16 
particles per cell for both electrons and ions, requiring a total of 5000 CPU hours on 
100 processors to simulate more than 1500 laser periods of the interaction. The code is 
fully parallelized and the simulations were performed at the CINECA supercomputing 
facility in Bologna (Italy). 

In the following, lengths are given in units of A, times in units of Tl = A/c = 27r/o>, 
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electric and magnetic fields in units of E = m e ujc/e, and densities in units of 
n c = m e u 2 /ATie 2 . For A = 1 /im, E = 3.213 x 10 10 V cm" 1 = 107.1 MGauss and 
n c = 1.11 x 10 21 cm -3 . The dimensionless parameter at,, giving the peak field amplitude 
of the laser pulse normalized to E , is related to the laser intensity / and wavelength 
by a L = 0.85(/A 2 /10 18 Wcni^/im 2 ) 1 / 2 . 

In all the 2D simulations reported below, the plasma is inhomogeneous along the x 
axis, i.e. in the direction of propagation of the laser pulse. The electron density profile 
rises linearly from zero value at x = 25A to the peak value no = 0.1n c at x = 425A, and 
then remains uniform. The pulse duration was either 150 or 300 Tl, corresponding 
to 0.5 and 1 ps, respectively, for A = 1 /im. 

The laser pulse was S-polarized, i.e. the electric field of the laser pulse was in 
the z direction perpendicular to the simulation plane. In the following we restrict 
the discussion to the S'-polarization case which has some advantages for the data 
analysis and visualization (for instance, the space-charge field generated in radial (y) 
direction during self-channeling is separated by the electromagnetic field E z which is 
representative of the pulse evolution). It is known, however, that at high intensity the 
details of nonlinear effects in pulse propagation depend on the polarization, leading to 
differences between the 5"- and P-polarization cases in 2D geometry and to asymmetry 
effects in 3D for what concerns self-focusing [23] and also to differences in the type 
and stability of solitons and vortices [21 El and references therein]. A preliminary 
simulation performed for P-polarization showed slight, but no substantial differences for 
what concerns the early self-channeling evolution which we discuss in section 13.11 The 
discussion of the effect of different polarizations on the coherent structures generation 
and evolution is more involved and will be addressed in future work. 

3. Results 

To illustrate the variety of nonlinear effects observed in the simulation results, Figure [U 
shows snapshots at t = 10 3 Tt of the ion density (rij) and the electric field of the laser 
pulse (E z ) over nearly the whole length of the plasma, for a simulation with ai = 2.7 and 
tl = 330Tt. Figure [1] contains most of the prominent features we observed throughout 
the set of our simulations, which may be summarized as follows. 

In the low-density region, the laser pulse bores a single charge-displacement channel, 
which in the higher density region breaks up into three main channels and a few 
of secondary, narrow filaments. In the following (see section 13.11) we trace back the 
appearance of the "trifurcated" channel to the effects of radial ion acceleration which 
lead to the "breaking" of the channel walls. 

Different types of electromagnetic structures are observed in regions of different 
density. In the lower density region (approximately between x = 100 and x = 150 in 
Figured]) a pattern of fields with approximate axial symmetry is observed. The detailed 
analysis of the electric and magnetic fields, including an estimate of their characteristic 
frequency from the simulation (see section 13.41) shows that this type of structures 
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combines both features of low-frequency electromagnetic post-solitons or "cavitons" 
and steady current vortices. In the higher density region a number of slowly-evolving 
field structures, either appearing as "solitary" structures or organized into patterns, are 
observed both outside the main low-density channels and inside the latter. There is 
some experimental indication of the growth of regularly spaced field structures into the 
main channel [23] . 

3.1. Ion and electric field dynamics following self- channeling 

For intensities up to ~ 2, in the early stage of the interaction the laser pulse bores a 
regular charge-displacement channel in the inhomogeneous region of the plasma, i.e. at 
densities n e < 0.1n c . This is the case for the simulation of Figure [2] (ai = 2, tl = 300Tl, 
transverse width r^ = 4 A) which shows a snapshot of the ion density rij and the electric 
field components E z and E y (results from this simulation are also reported in [22]). The 
laser pulse undergoes self-focusing as indicated both by the reduction of its transverse 
radius to ~ 3 A and the increase of its amplitude by a factor ~ 1.2. 

In the leading edge of the channel the transverse field E y is in the outward direction 
from the axis, indicating that the channel is positively charged due to the radial 
expulsion of electrons. In the trailing part of the pulse, the radial profile of E y changes 
qualitatively, as two ambipolar fronts appear on each side of the channel. On the inner 
side of the ambipolar fronts E y now points in the inward direction, i.e. towards the 
axis. The onset of an "inversion" in the radial field has been noticed in experimental 
investigations of channel dynamics [22] . 

3.2. One- dimensional modeling and the electric field "echo" 

The dynamics leading to the evolution of the radial electric field can be studied in detail 
using an one-dimensional, electrostatic PIC model where the laser pulse action is taken 
into account only via the ponderomotive force. The model assumes a non-evolving radial 
profile of the laser pulse and cylindrical symmetry, taking only the radial, cycle-averaged 
dynamics of electron and ions into account. Details about the model and its results are 
reported elsewhere [25]. Here we focus on the most prominent features of electric field 
dynamics. 

Figure [2] shows snapshots of the radial electric field E r and the ion density rij 
at various times, for a ID simulation in the same regime of 2D electromagnetic runs. 
Initially, the ponderomotive force F p pushes electron away from the axis, creating a 
back-holding space-charge field which is found to balance F p almost exactly. At the end 
of the pulse, when F p = 0, E r has almost vanished. However, E r appears back at a later 
time, with an ambipolar profile very similar to that observed in the 2D simulations. 
This "echo" effect originates from the ion dynamics of ions which are accelerated by 
the electric force ZeE r = ZF P during the laser pulse. The spatial profile of F p is such 
that the ions are focused towards a very narrow region at the edge of the channel, 
producing a very sharp spike of the ion density and leading to hydrodynamical breaking 
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as the fastest ions overturn the slowest ones. Looking at the profile of the ion density 
we observe that the latter may be said to "break" in literal meaning, as a secondary 
density spike moving outwards is formed. The process is also accompanied by strong 
heating of electrons near the breaking point, leading to the appearance of an ambipolar 
sheath field around the density spike. The negative field is strong enough to slow down 
and invert the velocity of the slowest ions, which are directed back to the axis where 
they are found to form a local density maximum at later times. 

3.3. Laser beam breakup 

A simple analytical model shows that the time required for the ions in the channel 
to reach the "breaking" point is proportional to the channel radius and inversely 
proportional to the laser field amplitude [25]. For high intensities, the "breaking" 
effect due to ion acceleration may occur early during the laser pulse, i.e. when the 
electromagnetic energy density inside the channel is very high, and cause a fast, strong 
variation of the density at the edge of the channel. In turn, this may affect the 
propagation of the laser pulse, similarly to what would happen in a wave guide where a 
sudden "leak" in its walls occurs. A possible signature of this effect is the appearance of 
two secondary beams, propagating in oblique direction, and originating near the point 
where the breaking of the channel walls occurs, as can be observed in Figure [3j 

From the "leaking waveguide" picture we roughly estimate these secondary beams 
to propagate at an angle 9 with respect to the axis given by tan 9 ~ k y /k x , where 
k y ~ Ti/d is the transverse wavevector of the guided mode, d is the local channel 



diameter, and k x ~ ^Juj 2 /c 2 — k 2 . In this estimate the pulse in channel is modelized 
as a TE mode of lowest order in a square guide. From the simulation result we get 
tan# ~ 0.065, while being d ~ 7A we obtain k y /k x ~ (n/7\)/(2n/\) = 0.071. 

3.4- Slowly-varying electromagnetic structures 

As already noticed in Figure [1] an impressive number of localized, slowly varying 
structures is generated in the interaction. In the denser plasma region, the several 
small-scale structures whose most evident signature is a strong depression in the plasma 
density are likely to be rather similar to the so-called post-solitons (T2J [14] , having zero 
propagation velocity and slowly expanding due to ion acceleration driven by the internal 
radiation pressure. They may be described as small cavities trapping electromagnetic 
radiation whose frequency is less than the plasma frequency of the surrounding plasma 
(hence they may be also appropriately named as electromagnetic "cavitons"). We notice 
that we do not observe a drift of such structures towards the low-density region. This 
difference from the observations of Ref. [11] might be ascribed to the smoother electron 
gradient in our case. 

The regular structures, forming an axially symmetrical row, observed in the low 
density region near the plasma boundary (far left side in Figure [1]) have indeed features 
which are similar both to electromagnetic cavitons and magnetic vortices. This "dual" 
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nature can be observed in Figure HI which shows the components of the fields E z and 
B z perpendicular to the simulation plane as a contour plot and the components in the 
(x, y) plane as a vector plot. By analyzing the frequency spectrum of the fields inside the 
density depression, we find that the fields E z , B x and B y are oscillating at a frequency of 
approximately O.Icj, lower than the local value of the plasma frequency (for unperturbed 
plasma) uj p ~ 0.15a;. Qualitatively, the oscillating fields are similar to those of the lowest 
TM resonant mode in a cylindrical cavity. 

The frequency analysis of E x , E y and B z shows that these field components are 
quasi-static, their spectrum being peaked around zero frequency. The electric field 
components E x and E y are in radial direction with respect to the axis of the structure, 
as it is expected for a cavity expanding under the action of the radiation pressure of 
the trapped radiation. The static magnetic field component B z is associated to current 
rings flowing around the axis of the structure. 

Apart from being associated to "post-soliton"-like structures, the fact that the 
magnetic vortices form a symmetrical row and are localized near the boundary of the 
channel make them different from those observed in the wake of a much shorter laser 
pulse, for which the creation of a low-density channel does not occur, and which seem 
to form an antisymmetrical row (TJ [26]. It is nevertheless possible that the current 
filamentation instability discussed in Ref. [26] plays a role in vortex formation also in 
the present case. In the early stage we observe a strong electron current in the main 
channel and two narrow return current sheets just outside the channel boundaries; later, 
the current layers seem to bend locally forming vortices around magnetic field maxima. 

The axial symmetry of these particular structures suggests that in "realistic" 3D 
geometry they may have a toroidal or "donut" shape. To get an impression of such 
a 3D structure one should imagine to rotate the field patterns of the 2D simulations 
around the x axis. This particular type of coherent structure would be characterized by 
azimuthal components of E (oscillating) and B (quasi-static) directed along the torus 
circumference, a solenoidal and oscillating magnetic field coiled up round the torus, 
and by an electrostatic fields component perpendicular to the torus surface. The 3D 
soliton discussed in Ref . [13] has a toroidal magnetic field and a poloidal magnetic field; 
however, in our case we have no clear indication of the charge oscillations inside the 
solitons observed in Ref. [13] . 

Inside the main and secondary low-density channels generated in the denser region 
of the plasma, the growth of field patterns which are less regular than those of Figure HI 
but qualitatively similar can be observed. Figure [5] give details of their evolution. We 
observe a tendency of this type of structures to grow inside the channels and to be 
correlated with rippling and bending of the channel walls. Theoretical work will be 
required to address the physics of formation of such structure patterns. 
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4. Conclusions 

The main results emerging from the series of 2D PIC simulations reported in the present 
paper may be summarized as follows. Ion acceleration due to the space-charge field in the 
channel drilled by the laser pulse leads to hydrodynamical breaking of the plasma profile 
at the channel walls. Two side effects of the ion-driven "breaking" have been identified: 
a change in the radial profile of the electrostatic field (including a sort of "echo" effect for 
pulses shorter than the breaking time) and a breakup of "long" laser pulses due to the 
sudden "leak" generated in the channel walls. The evolution of coherent, slowly- varying 
field structures has been monitored in time up to thousand of laser cycles, corresponding 
to several picoseconds in "real" experiments. Patterns of multi-peak structures appear 
inside low-density channels, and the formation of structures having both oscillating and 
static field components with an hybrid soliton-vortex nature has been observed. These 
results, and the perspective of experimental investigations of such field patterns, support 
the view of relativistic "laser plasmas" as environments showing an high degree of self- 
organization and a wealth of coherent structures, which are thus of great interest for 
the physics of nonlinear systems. 
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Figure 1. Contours of the electric field of the laser pulse E z and the ion density rii from 
a 2D PIC simulation, at t = 10007^, showing the self-channeling and filamcntation of 
the laser pulse and the generation of isolated solitary structures and of field patterns. 
The laser pulse propagates from left to right. In the E z frame, the contour levels in 
the leftmost region (50 < x < 250) has been rescaled by a factor of 3 to show the 
presence of field structures with relatively low amplitudes. The laser pulse parameters 
are ql — 2.7, = 8A, and tl = 3307l. 



Figure 2. Simulations results addressing electric field dynamics following self- 
channeling, showing the transition in the radial field profile [22]. Left column: 2D 
PIC results. Top frame: ion density (rii) and electric field components (E z and E y ) 
at t = 6002l. Bottom frame: lineout of E y (blue) and rii (red) along the y-axis at 
two different re-positions. Parameters are = 2, tl = 3001l, tl — 4 A. Right 
column: snapshots at various times of radial electric field E r (blue, thick line) and ion 
density n; (red, dash-dotted line), and the phase space distributions of ions fi(r,p r ) 
and electrons / e (f, p r ) from ID simulations using a ponderomotive, electrostatic model 
[2"5] , Parameters are = 2.7, n e /n c = 0.01, rj, = 7.5A, tx = 330T^. 



Figure 3. Evolution of the laser field E z at different times showing the breakup of the 
laser pulse into three main beams. The two secondary beams propagating in oblique 
direction originate from near the location of the "breaking" of the channel walls. The 
laser pulse parameters are = 2.7, Tl = 8A, and tl = 150Tz,. 



Figure 4. (Anti-)symmctrical row of slowly-varying structures in the low density 
region of the plasma at t = 6257^. The left column shows the fields E z (contour plot) 
and B x + B y (vector plot) oscillating at a frequency ~ O.lw. The right column shows 
the quasi-static fields B z (contour plot) and E x + E y (vector plot). The laser pulse 
parameters are = 2.7, = 8A, and tl — 1507l. 



Figure 5. Detail of the evolution of field structures in the denser region of the plasma, 
from the simulation of Figure [TJ The laser pulse parameters are = 2.7, rt = 8A, 
and t l = 300T L . 
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